###Bellman equation for working periods
function Update_valfunc(prim::Primitives, est::Estimands, res::Results, i_t::Int64, i_θ::Int64, i_e::Int64, utilities::Vector{Matrix{Float64}}; rec::Int64 = 0, decomp::Int64 = 0)
    @unpack r, η, β, cfloor, θ_grid, nθ, a_grid, na, t_grid, nt, e_grid, ne, x_grid, nx, r_d, d_grid, nd, prob_forb, p_min, p_max = prim
    @unpack γ_grid, λ_0, λ_1, λ_2, λ_3, λ_4, ξ, σ_l, π_0, π_1, π_2, π_3, π_4, σ_ζ, ψ_0, ψ_1, ψ_2 = est

    t = t_grid[i_t] #age

    if i_t == nt #in terminal period
        for i_x = 1:nx, i_a = 1:na, i_d = 1:nd  #loop over state space
            θ, x, a, d = θ_grid[i_θ], x_grid[i_x], a_grid[i_a], d_grid[i_d] #state variables
            
            #leisure utility and income
            util_leisure = λ_0 + λ_1*θ + λ_2*t + λ_3*(t-50)*(t>50) + λ_4 * t^3 
            incwage = exp(γ_grid[i_e, 1] + γ_grid[i_e, 2] * θ + γ_grid[i_e, 3] * x + γ_grid[i_e, 4]*x^2)

            #standard payment if still debt left
            pay_left = 1
            pay_stand = 0.0
            if pay_left>0
                pay_stand = payment(d, r_d, pay_left)

                #bounding
                pay_stand = max(pay_stand, p_min)
                pay_stand = min(pay_stand, p_max)
                pay_stand = min(pay_stand, d * (1+r_d)) #bound above by debt remaining
            end
            
            #a' must be zero, so just solve for consumption by eating everything. Including forebearance state.
            c_emp = incwage - tax(incwage) + (1+r)*a - pay_stand
            c_unemp = (1+r)*a - pay_stand 
            c_unemp_forb = (1+r)*a 

            #update value functions 
            res.val_func_emp[i_θ, i_e, i_x, i_a, i_t, i_d] = utility(max(c_emp, cfloor), η) 
            res.val_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = utility(max(c_unemp, cfloor), η) 
            res.val_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = utility(max(c_unemp_forb, cfloor), η) 
            res.val_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = utility(max(c_unemp, cfloor), η) + util_leisure
            res.val_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = utility(max(c_unemp_forb, cfloor), η) + util_leisure

            #update consumption policy functions
            res.pol_func_emp[i_θ, i_e, i_x, i_a, i_t, i_d] = c_emp
            res.pol_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = c_unemp
            res.pol_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = c_unemp_forb
            res.pol_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = c_unemp
            res.pol_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = c_unemp_forb
        end
    elseif i_t!=nt #not in terminal period!
        for i_x = 1:nx, i_a = 1:na, i_d = 1:nd  #loop over state space
            θ, x, a, d = θ_grid[i_θ], x_grid[i_x], a_grid[i_a], d_grid[i_d] #state variables
            x_check = i_x + 2*(i_e == 2) + 4 * (i_e == 3)

            #skip if experience grid is higher than age grid, since this is impossible. 
            if x_check > i_t
                continue
            end

            #assume loans paid off by age 42. If not, you're in trouble!!
            if i_t>=34 && i_d!=1
                continue
            end

            ###figure out loan payment###
            
            #approximate number of payment periods left. Start at 19 if some college
            periods_left = 10 - (i_t - 1) + 1
            periods_left += 1*(i_e==2) #update if somecoll
            periods_left += 3*(i_e==3) #update if college
            
            #bound
            periods_left = max(periods_left, 1)
            periods_left = min(periods_left, 10)
            
            #compute standard payment and bound
            pay_stand = payment(d, r_d, periods_left)
            pay_stand = max(pay_stand, p_min)
            pay_stand = min(pay_stand, p_max)
            pay_stand = min(pay_stand, d * (1+r_d)) #bound above by debt remaining

            ###figure out bounds for consumption choices###
            incwage = exp(γ_grid[i_e, 1] + γ_grid[i_e, 2] * θ + γ_grid[i_e, 3] * x + γ_grid[i_e, 4]*x^2)
    
            #minimum: can't consume less than the floor!
            cmin = cfloor
            
            #max consumption levels: bounded by a'>=-cfloor
            cmax_emp = incwage - tax(incwage) + (1+r)*a - pay_stand + cfloor 
            cmax_unemp = (1+r)*a - pay_stand + cfloor
            cmax_unemp_forb = (1+r)*a + cfloor

            ####start solving policy functions and updating value functions#####

            #define criterion functions
            f_emp(c) = Optimal_Consumption(c, prim, est, res, i_θ, i_e, i_x, i_a, i_t, i_d, 2, utilities; rec=rec, decomp=decomp) * -1
            f_unemp(c) = Optimal_Consumption(c, prim, est, res, i_θ, i_e, i_x, i_a, i_t, i_d, 1, utilities; rec=rec, decomp=decomp) * -1
            f_unemp_forb(c) = Optimal_Consumption(c, prim, est, res, i_θ, i_e, i_x, i_a, i_t, i_d, 1, utilities; rec=rec, decomp=decomp, forb=1) * -1
            f_nilf(c) = Optimal_Consumption(c, prim, est, res, i_θ, i_e, i_x, i_a, i_t, i_d, 0, utilities; rec=rec, decomp=decomp) * -1
            f_nilf_forb(c) = Optimal_Consumption(c, prim, est, res, i_θ, i_e, i_x, i_a, i_t, i_d, 0, utilities; rec=rec, decomp=decomp, forb=1) * -1
            
            #employed
            if cmax_emp<cfloor #borrowing-constrained -- gov't steps in and provides cfloor. Agent pays off as much debt as possible.  
                #update value/policy function
                res.val_func_emp[i_θ, i_e, i_x, i_a, i_t, i_d] = f_emp(0.0)
                res.pol_func_emp[i_θ, i_e, i_x, i_a, i_t, i_d] = 0.0
            elseif cmax_emp>cfloor #not borrowing-constrained -- now we can choose stuff!
                #update value/policy function
                opt_c_emp = Optim.optimize(f_emp, cmin, cmax_emp; abs_tol = 1e-4, rel_tol = 1e-2)
                res.val_func_emp[i_θ, i_e, i_x, i_a, i_t, i_d] = opt_c_emp.minimum * -1
                res.pol_func_emp[i_θ, i_e, i_x, i_a, i_t, i_d] = opt_c_emp.minimizer
            end

            #unemployed/NILF, not in forebearance
            if cmax_unemp<cfloor #borrowing-constrained -- gov't steps in and provides cfloor. Agent pays off as much debt as possible.
                #update value/policy function for unemp
                res.val_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = f_unemp(0.0)
                res.pol_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = 0.0

                #update value/policy function for nilf
                res.val_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = f_nilf(0.0)               
                res.pol_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = 0.0
            elseif cmax_unemp>cfloor #not borrowing-constrained                
                #update value/policy function for unemp
                opt_c_unemp = Optim.optimize(f_unemp, cmin, cmax_unemp; abs_tol = 1e-4, rel_tol = 1e-2)
                res.val_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = opt_c_unemp.minimum * -1
                res.pol_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = opt_c_unemp.minimizer
                
                #update value/policy function for nilf
                opt_c_nilf = Optim.optimize(f_nilf, cmin, cmax_unemp; abs_tol = 1e-4, rel_tol = 1e-2)
                res.val_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = opt_c_nilf.minimum * -1                
                res.pol_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 1] = opt_c_nilf.minimizer
            end

            #unemployed in forebearance
            if i_d == 1 || i_t>=34 #irrelevant if no debt or too old
                res.val_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = res.val_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 1]
                res.pol_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = res.pol_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 1]

                #update value/policy function for nilf
                res.val_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = res.val_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 1] 
                res.pol_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = res.pol_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 1] 
            elseif cmax_unemp_forb<cfloor #borrowing-constrained -- gov't steps in and provides cfloor. Agent pays off as much debt as possible.
                #update value/policy function for unemp
                res.val_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = f_unemp_forb(0.0)
                res.pol_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = 0.0

                #update value/policy function for nilf
                res.val_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = f_nilf_forb(0.0)               
                res.pol_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = 0.0
            elseif cmax_unemp_forb>cfloor #not borrowing-constrained     
                #update value/policy function for unemp
                opt_c_unemp_forb = Optim.optimize(f_unemp_forb, cmin, cmax_unemp_forb; abs_tol = 1e-4, rel_tol = 1e-2)
                res.val_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = opt_c_unemp_forb.minimum * -1
                res.pol_func_unemp[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = opt_c_unemp_forb.minimizer
                
                #update value/policy function for nilf
                opt_c_nilf_forb = Optim.optimize(f_nilf_forb, cmin, cmax_unemp_forb; abs_tol = 1e-4, rel_tol = 1e-2)
                res.val_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = opt_c_nilf_forb.minimum * -1                
                res.pol_func_nilf[i_θ, i_e, i_x, i_a, i_t, i_d, 2] = opt_c_nilf_forb.minimizer
            end
        end
    end
end


function Optimal_Consumption(c::Float64, prim::Primitives, est::Estimands, res::Results, i_θ::Int64, i_e::Int64, i_x::Int64, i_a::Int64, i_t::Int64, i_d::Int64,
    l::Int64, utilities::Vector{Matrix{Float64}}; rec::Int64 = 0, decomp::Int64 = 0, forb::Int64 = 0)
    @unpack r, η, β, cfloor, θ_grid, nθ, a_grid, na, t_grid, nt, e_grid, ne, x_grid, nx, r_d, d_grid, nd, prob_forb, p_min, p_max = prim
    @unpack γ_grid, λ_0, λ_1, λ_2, λ_3, λ_4, ξ, σ_l, π_0, π_1, π_2, π_3, π_4, π_5, π_6 = est

    #unpack state space and relevant variables
    θ, x, a, d = θ_grid[i_θ], x_grid[i_x], a_grid[i_a], d_grid[i_d] #state variables
    t = t_grid[i_t] #age
    incwage = exp(γ_grid[i_e, 1] + γ_grid[i_e, 2] * θ + γ_grid[i_e, 3] * x + γ_grid[i_e, 4]*x^2)
    util_leisure = λ_0 + λ_1*θ + λ_2*t + λ_3*(t-50)*(t>50) + λ_4 * t^3
    
    #compute loan payment
    periods_left = 10 - (i_t - 1) + 1
    periods_left += 1*(i_e==2) #update if somecoll
    periods_left += 3*(i_e==3) #update if college
    
    #bound
    periods_left = max(periods_left, 1)
    periods_left = min(periods_left, 10)
    
    #compute standard payment and bound
    pay_stand = payment(d, r_d, periods_left)
    pay_stand = max(pay_stand, p_min) #bound below
    pay_stand = min(pay_stand, p_max) #bound above by max payment
    pay_stand = min(pay_stand, d * (1+r_d)) #bound above by debt remaining

    ###construct flow utility###
    util_flow = utility(max(c, cfloor), η) + util_leisure*(l==0) #easy!

    ###construct continuation payoff###
    ###first: figure out next-period indices
    i_tp = i_t + 1 #one year older
    i_xp = i_x + 1*(l==2) #one more year exp if worked
    i_θp, i_ep = i_θ, i_e #these don't change
    xp, ep = x_grid[i_xp], e_grid[i_ep]

    #next-period assets, student debt, and payment status
    ap_val, dp_val = 0.0, 0.0
    i_ap, i_dp = 0, 0

    if c == 0.0 #borrowing constrained -- here things are interesting
        ap_val = -cfloor #definitional -- borrowing constraint binds. More interesting is what happens to student loans
        dp_val = d * (1+r_d) - pay_stand*(1-forb) #evolution of debt -- appreciates if in forebearance

        #delinquincy check. This happens if can't pay off loans in full, even with c = cfloor and a' = -cfloor, and not in forebearance
        if (1+r)*a - pay_stand<0 && forb == 0
            pay_del = max((1+r)*a, 0)
            dp_val = (d*(1 + r_d) - pay_del) + 0.185 * (pay_stand - pay_del)
        end
        i_dp = get_index(dp_val, d_grid)
        i_ap = get_index(ap_val, a_grid)
    else #not borrowing constrained -- everything proceeds as usual
        ap_val = (incwage - tax(incwage)) * (l==2) + (1+r)*a - pay_stand*(1-forb) - c  #asset law of motion
        dp_val = d * (1+r_d) - pay_stand*(1-forb) #evolution of debt -- appreciates if in forebearance
        i_ap = get_index(ap_val, a_grid)
        i_dp = get_index(dp_val, d_grid)
    end

    ###construct threshold for LFP shocks
    #interpolated functions
    interp_val_emp = interpolate(res.val_func_emp[i_θp, i_ep, i_xp, :, i_tp, :], BSpline(Linear()))
    interp_val_unemp = interpolate(res.val_func_unemp[i_θp, i_ep, i_xp, :, i_tp, :, 1], BSpline(Linear()))
    interp_val_unemp_forb = interpolate(res.val_func_unemp[i_θp, i_ep, i_xp, :, i_tp, :, 2], BSpline(Linear()))
    interp_val_nilf = interpolate(res.val_func_nilf[i_θp, i_ep, i_xp, :, i_tp, :, 1], BSpline(Linear()))
    interp_val_nilf_forb = interpolate(res.val_func_nilf[i_θp, i_ep, i_xp, :, i_tp, :, 2], BSpline(Linear()))

    #probability of employment next period
    #adjust employment probability if in recession
    rec_emp_shocks = utilities[7]
    π_0_r = π_0
    π_2_r = π_2
    π_3_r = π_3
    if rec==1 && i_t<6 && (decomp == 0 || decomp == 1)
        π_0_r = π_0 * rec_emp_shocks[i_t, 1]
        π_2_r = π_2 * rec_emp_shocks[i_t, 2]
        π_3_r = π_3 * rec_emp_shocks[i_t, 3]
    end
    prob_emp = 𝚽(π_0_r + π_1 * (l==2) + π_2_r * (ep==2) + π_3_r * (ep==3) + π_4 * xp + π_5 * xp * (ep==2) + π_6 * xp * (ep==2))

    #values of working, not woring
    val_work = prob_emp * interp_val_emp(i_ap, i_dp) + (1-prob_emp) * (interp_val_unemp(i_ap, i_dp) * (1-prob_forb) +  interp_val_unemp_forb(i_ap, i_dp) * prob_forb)
    val_nwork = interp_val_nilf(i_ap, i_dp) * (1-prob_forb) + interp_val_nilf_forb(i_ap, i_dp) * prob_forb
    
    #compute threshold
    ε_threshold = val_nwork - val_work - ξ
    
    #if didn't work in current period, now has to pay the switching cost to work instead of having it taken off
    if l == 0
        ε_threshold += (2*ξ)
    end

    ###assmelbe expected max utility applying Mills ratio
    util_cont = cdf(Normal(0, σ_l), ε_threshold) * val_nwork + (1-cdf(Normal(0, σ_l), ε_threshold)) * (val_work + Mills(ε_threshold, 0.0, σ_l))
    val = util_flow + β * util_cont
    val
end

####Parent Bellman equation
function Valfunc_parent(prim::Primitives, est::Estimands, res::Results, utilities::Vector{Matrix{Float64}}; 
    rec::Int64 = 0, cfact::Int64 = 0, pbr::Int64 = 1, dlim::Int64 = 0, decomp::Int64 = 0, rec_pol::Int64 = 0)
    @unpack r, η, β, cfloor, dbar, I_grid, nI, W_grid, nW, θ_grid, nθ, a_grid, na, t_grid, nt, e_grid, ne, x_grid, nx, r_d, d_grid, nd, prob_forb, thresh_sub, dbar_baseline, dbar_nolim, dbar_rec = prim
    @unpack γ_grid, λ_0, λ_1, λ_2, λ_3, λ_4, ξ, σ_l, π_0, π_1, π_2, π_3, π_4, π_5, ω, α, ψ_0, ψ_1, ψ_2, σ_ζ, ψ_3, ψ_4 = est

    #yank things out of utilities package
    educ_costs, grad_probs, grants_func, grants_merit_func = utilities[1], utilities[2], utilities[3], utilities[4]
    d_opts = utilities[9]

    #setup
    γ_eul = Base.MathConstants.eulergamma

    #debt forgiveness
    d_reduce = 0 #debt forgiveness

    #recession: Stafford expansion, unless in alternate recession with no educational policy changes
    if rec_pol == 1 && pbr == 1 && dlim==0 && (decomp==0 || decomp == 3)
        dbar = dbar_rec
    end

    #hard-coding of borrowing constraints
    if dlim != 0 && pbr == 1
        dbar = d_opts[dlim, 1]        
    end

    @sync @distributed for i_θ = 1:nθ #loop over ability
        θ = θ_grid[i_θ] #level value

        #interpolated value functions, employed
        interp_v_hs_emp = interpolate(res.val_func_emp[i_θ, 1, 1, :, 1, :], BSpline(Linear())) #hs
        interp_v_drop_emp = interpolate(res.val_func_emp[i_θ, 1, 1, :, 2, :], BSpline(Linear())) #dropout
        interp_v_sc_emp = interpolate(res.val_func_emp[i_θ, 2, 1, :, 3, :], BSpline(Linear())) #sc
        interp_v_coll_emp = interpolate(res.val_func_emp[i_θ, 3, 1, :, 5, :], BSpline(Linear())) #college

        #interpolated value functions, unemployed
        interp_v_hs_unemp = interpolate(res.val_func_unemp[i_θ, 1, 1, :, 1, :, 1], BSpline(Linear())) #hs
        interp_v_drop_unemp = interpolate(res.val_func_unemp[i_θ, 1, 1, :, 2, :, 1], BSpline(Linear())) #dropout
        interp_v_sc_unemp = interpolate(res.val_func_unemp[i_θ, 2, 1, :, 3, :, 1], BSpline(Linear())) #sc
        interp_v_coll_unemp = interpolate(res.val_func_unemp[i_θ, 3, 1, :, 5, :, 1], BSpline(Linear())) #college
        
        #interpolated value functions, nilf
        interp_v_hs_nilf = interpolate(res.val_func_nilf[i_θ, 1, 1, :, 1, :, 1], BSpline(Linear())) #hs
        interp_v_drop_nilf = interpolate(res.val_func_nilf[i_θ, 1, 1, :, 2, :, 1], BSpline(Linear())) #dropout
        interp_v_sc_nilf = interpolate(res.val_func_nilf[i_θ, 2, 1, :, 3, :, 1], BSpline(Linear())) #sc
        interp_v_coll_nilf = interpolate(res.val_func_nilf[i_θ, 3, 1, :, 5, :, 1], BSpline(Linear())) #college

        ###same thing for forebearance state
        interp_v_drop_unemp_forb = interpolate(res.val_func_unemp[i_θ, 1, 1, :, 2, :, 2], BSpline(Linear())) #dropout
        interp_v_sc_unemp_forb = interpolate(res.val_func_unemp[i_θ, 2, 1, :, 3, :, 2], BSpline(Linear())) #sc
        interp_v_coll_unemp_forb = interpolate(res.val_func_unemp[i_θ, 3, 1, :, 5, :, 2], BSpline(Linear())) #college
        
        interp_v_drop_nilf_forb = interpolate(res.val_func_nilf[i_θ, 1, 1, :, 2, :, 2], BSpline(Linear())) #dropout
        interp_v_sc_nilf_forb = interpolate(res.val_func_nilf[i_θ, 2, 1, :, 3, :, 2], BSpline(Linear())) #sc
        interp_v_coll_nilf_forb = interpolate(res.val_func_nilf[i_θ, 3, 1, :, 5, :, 2], BSpline(Linear())) #college

        #begin loop over parent resources
        for i_I = 1:nI, i_W = 1:nW 
            educ_costs_temp = copy(educ_costs)
            I, W = I_grid[i_I], W_grid[i_W] #state variables
            candidate_τI, candidate_τW, candidate_τ = 0, 0, 0 #candidate transfers
            candidate_max = -1e10 #new candiate max

            #modify tuition prices in recessions
            if rec == 1 && pbr == 1 && (decomp==0)
                educ_costs_temp[2] *= 1.00
                educ_costs_temp[3] *= 1.08
                educ_costs_temp[4] *= 1.03
                educ_costs_temp[5] *= 1.03
            end

            #finanical aid
            grants_aid = zeros(5)
            grants_fed = zeros(5)

            for s = 1:4
                grants_aid[s+1] = Grants(I, s, grants_func, W)[1]
                grants_fed[s+1] = Grants(I, s, grants_func, W)[2]

                if cfact>0 && pbr==1 #counterfactual Pell extension
                    grants_aid[s+1] = Grants(I, s, grants_func, W; cfact = cfact)[1]
                    grants_fed[s+1] = Grants(I, s, grants_func, W; cfact = cfact)[2]
                end
            end

            #merit grants
            grants_merit = zeros(5)
            for s = 1:4
                grants_merit[s+1] = Grants_merit(θ, s, grants_merit_func)
            end

            #uupdate pell receipts in recession
            if I<0.820 && rec_pol == 1 && pbr==1 && (decomp==0 || decomp == 3)
                #updating policy
                for s = 1:4
                    temp = grants_aid[s+1] - grants_fed[s+1] + (grants_fed[s+1]*1.25)
                    temp2 = grants_fed[s+1] * 1.25
                    grants_aid[s+1] = temp
                    grants_fed[s+1] = temp2
                end
            end

            #total grants
            grants_total = grants_aid .+ grants_merit

            #####loop over potential transfer combinations#####
            τI_grid, τW_grid = collect(0:0.025:I), collect(0:0.025:W) #2.5k increments
            nτI, nτW = length(τI_grid), length(τW_grid)

            for i_τI = 1:nτI, i_τW = 1:nτW #loop over transfers
                τI, τW = τI_grid[i_τI], τW_grid[i_τW] #transfer breakdown
                τ = τI + τW #trasnsfer offered to kid

                if τ>4.0 #upper bound on total transfer
                    continue
                end

                ###Start computing kid utility for each schooling option
                ns = length(educ_costs_temp) #number of schooling options
                val_kid = zeros(ns)

                #expected value of hiugh school option, given assets and debt
                a_start_hs = 0.0 #you're on your own!
                d_start_hs = 0.0 #no debt, though!
                a_hs_ind = get_index(a_start_hs, a_grid)
                d_hs_ind = get_index(d_start_hs, d_grid)

                prob_emp_hs = 𝚽(π_0)
                prob_emp_drop = 𝚽(π_0)
                prob_emp_sc = 𝚽(π_0 + π_2)
                prob_emp_coll = 𝚽(π_0 + π_3)
                rec_emp_shocks = utilities[7]            

                #adjust initial employment probabilities if in recession and parents paying attention
                if rec == 1 && pbr == 1 && (decomp == 0 || decomp == 1)
                    #prob_emp_hs = 𝚽(π_0 * rec_emp_shocks[1, 1])
                    #prob_emp_drop = 𝚽(π_0 * rec_emp_shocks[2, 1])
                    prob_emp_hs = 0.7367
                    prob_emp_drop = 0.7491
                    prob_emp_sc = 𝚽(π_0 * rec_emp_shocks[3, 1] + π_2 * rec_emp_shocks[3, 2])
                    prob_emp_coll = 𝚽(π_0 * rec_emp_shocks[5, 1] + π_3 * rec_emp_shocks[5, 3])
                end

                #compute expected HS utility. No need to worry about forebearance here, since there's no student debt
                val_work = prob_emp_hs * interp_v_hs_emp(a_hs_ind, d_hs_ind) + (1-prob_emp_hs) * interp_v_hs_unemp(a_hs_ind, d_hs_ind)
                val_nwork = interp_v_hs_nilf(a_hs_ind, d_hs_ind)

                #working threshold rule (no switching costs!) and expected max utility
                ε_l_thresh = val_nwork - val_work
                val_hs = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                val_kid[1] = val_hs

                ###loop over other schooling options -- here's where things get more complicated
                for i_s = 2:ns 
                    prob_drop, prob_sc, prob_coll = Grad_prob(θ, i_s-1, grad_probs) #probability of graduation outcomes

                    ###starting assets/debt for each possible educational outcome. Bounded by zero in either case
                    #dropout
                    a_start_drop = max((τ - (educ_costs_temp[i_s] - grants_total[i_s]))/4, 0)
                    d_start_drop = min((τ - (educ_costs_temp[i_s] - grants_total[i_s]))/4, 0) * -1 

                    #some college
                    a_start_sc = max((τ - (educ_costs_temp[i_s] - grants_total[i_s]))/2, 0)
                    d_start_sc = min((τ - (educ_costs_temp[i_s] - grants_total[i_s]))/2, 0) * -1 

                    #college
                    a_start_coll = max((τ - (educ_costs_temp[i_s] - grants_total[i_s])), 0)
                    d_start_coll = min((τ - (educ_costs_temp[i_s] - grants_total[i_s])), 0) * -1 

                    #sub/unsubsidized loan check. If parent income is too high, interest accumulates while in college.   
                    if I>thresh_sub
                        d_start_drop *=(1+r_d)
                        d_start_sc *= (1+r_d)^2
                        d_start_coll *= (1+r_d)^4
                    end

                    a_d_i = get_index(a_start_drop, a_grid)
                    a_s_i = get_index(a_start_sc, a_grid)
                    a_c_i = get_index(a_start_coll, a_grid)
                    d_d_i = get_index(d_start_drop, d_grid)
                    d_s_i = get_index(d_start_sc, d_grid)
                    d_c_i = get_index(d_start_coll, d_grid)

                    ###utility from dropout: threshold, then compute
                    val_work = prob_emp_drop * interp_v_drop_emp(a_d_i, d_d_i) + (1-prob_emp_drop) * (interp_v_drop_unemp(a_d_i, d_d_i) * (1-prob_forb) + interp_v_drop_unemp_forb(a_d_i, d_d_i) * prob_forb)
                    val_nwork = interp_v_drop_nilf(a_d_i, d_d_i) * (1-prob_forb) + interp_v_drop_nilf_forb(a_d_i, d_d_i) * prob_forb
                    ε_l_thresh = val_nwork - val_work
                    val_drop = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                    val_kid[i_s] += β * val_drop * prob_drop

                    #utility from some college
                    val_work = prob_emp_sc * interp_v_sc_emp(a_s_i, d_s_i) + (1-prob_emp_sc) * (interp_v_sc_unemp(a_s_i, d_s_i)*(1-prob_forb) + interp_v_sc_unemp_forb(a_s_i, d_s_i) * prob_forb)
                    val_nwork = interp_v_sc_nilf(a_s_i, d_s_i) * (1-prob_forb) + interp_v_sc_nilf_forb(a_s_i, d_s_i) * prob_forb
                    ε_l_thresh = val_nwork - val_work
                    val_sc = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                    val_kid[i_s] += β^2 * val_sc * prob_sc

                    #utilty from college
                    val_work = prob_emp_coll * interp_v_coll_emp(a_c_i, d_c_i) + (1-prob_emp_coll) * (interp_v_coll_unemp(a_c_i, d_c_i) * (1-prob_forb) + interp_v_coll_unemp_forb(a_c_i, d_c_i) * prob_forb)
                    val_nwork = interp_v_coll_nilf(a_c_i, d_c_i) * (1-prob_forb) + interp_v_coll_nilf_forb(a_c_i, d_c_i) * prob_forb
                    ε_l_thresh = val_nwork - val_work
                    val_coll = cdf(Normal(0, σ_l), ε_l_thresh) * val_nwork + (1-cdf(Normal(0, σ_l), ε_l_thresh)) * (val_work + Mills(ε_l_thresh, 0.0, σ_l))
                    val_kid[i_s] += β^4 * val_coll * prob_coll

                    #psychic cost of enrollment
                    d_exp = prob_drop * d_start_drop + prob_sc * d_start_sc + prob_coll * d_start_coll    
                    ψ_total = ψ_0 + ψ_1 * θ + ψ_2 * θ * (i_s==3 || i_s==5) + ψ_3 * (d_exp>0) + ψ_4 * d_exp
                    val_kid[i_s] += ψ_total

                    #check that we don't go over the borrowing constraint
                    debt_raw_max = τ - (educ_costs_temp[i_s] - grants_total[i_s])
                    if (debt_raw_max + d_reduce)<(-dbar) #over borrowing limit. Can't go over limit with forgiven debt
                        val_kid[i_s] = -Inf #not allowed to exceed borrowing constraint!
                    end                   
                end

                #expected kid utillity conditional on trasnfer
                e_val_kid = γ_eul * σ_ζ + σ_ζ * log(sum(exp.((1/σ_ζ).*val_kid)))

                ####construct parent utility
                ###initialize parent utiltiy based on altruistic benefit
                val_parent = α * e_val_kid 
                
                ###now add parent utility from resources
                val_kid_exp = exp.((1/σ_ζ).*val_kid)

                #high school -- parent keeps (1-f) τ of stuff
                prob_hs = val_kid_exp[1] / sum(val_kid_exp)
                val_parent += prob_hs * (log(I) + ω * log(W))

                #loop over other schooling options
                for i_s = 2:ns
                    
                    #probabiliiyt of selection
                    prob_s = val_kid_exp[i_s]/sum(val_kid_exp)
                    prob_drop, prob_sc, prob_coll = Grad_prob(θ, i_s-1, grad_probs) #probability of graduation outcomes

                    #parent utility over possible schooling outcomes
                    val_parent += prob_s * prob_drop * (log(I - τI/4) + ω * log(W - τW/4)) #if dropout, parent keeps some transfer
                    val_parent += prob_s * prob_sc * (log(I - τI/2) + ω * log(W - τW/2)) #otherwise lose it
                    val_parent += prob_s * prob_coll * (log(I - τI) + ω * log(W - τW)) #otherwise lose it
                end

                #check for new maximizer
                if val_parent>candidate_max
                    candidate_max = val_parent
                    candidate_τI, candidate_τW = τI, τW
                    candidate_τ = τ
                end
            end
            res.pol_func_p_I[i_I, i_W, i_θ] = candidate_τI #update policy function
            res.pol_func_p_W[i_I, i_W, i_θ] = candidate_τW #update policy function
            res.pol_func_p_τ[i_I, i_W, i_θ] = candidate_τ #update policy function
        end
    end
end
